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ABSTRACT 


It  is  well  known  [1],[2],(3]  that  if  T:Rn  -*■  Rn  is  a  Lipschitz 

continuous,  strongly  monotone  operator  and  X  is  a  closed  convex  set, 

* 

then  a  solution  x  €X  of  the  variational  inequality 
(x-x*)'f(x*)  >_  0,  V  x£X 

can  be  found  Iteratively  by  means  of  the  projection  method 

*k+i  “  V\  -  •  V*  • 

provided  the  stepsize  a  is  sufficiently  small.  We  show  that  the  same  is 
true  if  T  is  of  the  form 
T  -  A'TA, 

where  A:Rn  -*■  Rm  is  a  linear  mapping,  provided  T:Rm  Rm  is  Lipschitz 
continuous  and  strongly  monotone,  and  the  set  X  is  polyhedral.  This 
fact  is  used  to  construct  an  effective  algorithm  for  finding  a  network 
flow  which  satisfies  given  demand  constraints  and  is  positive  only  on 
paths  of  minimum  marginal  delay  or  travel  time. 
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Introduction 

•  • 

* 

We  consider  the  problem  of  finding  x  eX  satisfying  the  variational 
inequality 

(1)  (x-x*)'A'T(Ax*)  >_  0  ,  V  x£X 

where  X  is  a  nonempty  polyhedral  subset  of  Rn,  A  is  a  given  mxn  matrix, 
and  T:Rm  **■  Rm  is  a  nonlinear  operator  which  is  Lipschitz  continuous  and 
strongly  monotone  in  the  sense  that  there  exist  positive  scalars  L  and  X 
such  that  for  all  y^,  y^  in  the  set  Y  »•  {y|y«Ax,  x£X)  we  have 

C2)  |T(7l)  -  T(y2)|  <  L|yi  -  y2i 

(3)  [T(7l)  -  T(y2)]'  (yx  -  y2)  1  X|yx  -  y2|2  . 

In  the  relations  above  and  throughout  the  paper  all  vectors  are  considered 
to  be  column  vectors,  and  a  prime  denotes  transposition.  The  standard 
norm  in  Rn  is  denoted  | • | ,  i.e.  |x|  *  (x'x)1^2  for  all  xeRn. 

We  are  interested  in  the  projection  algorithm 

(A)  x^  -  P|[xk  -  aS-1  A’T(Ax^)] ,  xQeX 

where  a  >  0  is  a  stepsize  parameter,  S  is  a  fixed  positive  definite  symmetric 
matrix,  and,  for  any  zeRn,  P^(*)  denotes  the  unique  projection  of  z  on  the 
set  X  with  respect  to  the  norm  |  |  *  |  [  corresponding  to  S 

(5)  | |w| |  -  <w'Sw)1/2  ,  V  weRn. 

The  variational  inequality  (1)  arises  from  the  variational  inequality 

(6)  (y  -  y*)'  T(y*)  >  0  ,  V  yeY 

through  the  transformation 

(7)  y  -  Ax,  Y  -  AX  •  (y'y  -  Ax,  xeX} 
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It  Is  possible  to  employ  the  projection  algorithm 

(8)  ylc+l  "  PY^k  "  aQ’lT(yk)3  ’  yo£Y 

for  solving  (6)  where  a  is  a  positive  stepsize  parameter  and  Q  is  positive 
definite  symmetric.  It  has  been  shown  by  Sibony  [1]  that  if  T  is  Lipschitz 
continuous  and  strongly  monotone,  and  a  is  chosen  sufficiently  small  then 
the  sequence  {y^}  generated  by  (8)  converges  to  the  unique  solution  y  of 
(6) .  The  rate  of  convergence  is  typically  linear  although  a  superlinear 
convergence  rate  is  possible  in  exceptional  cases  [13].  Strong  monotonicity 
of  the  mapping  T  is  an  essential  assumption  for  these  results  to  hold.  Our 
motivation  for  considering  iteration  (4)  stems  from  the  fact  that  in  some 
cases  projection  on  the  set  Y  is  very  difficult  computationally  while 
projection  on  the  set  X  through  the  transformation  y  ■=  Ax  [c.f.  (7)]  may  be 
very  easy.  Under  these  circumstances  if  all  other  factors  are  equal,  the 
projection  method  (4)  is  much  more  efficient  than  the  method  (8) .  This 
situation  occurs  for  example  in  the  application  discussed  in  Section  3. 

A  potential  difficulty  with  the  transformation  idea  described 
above  is  that  the  mapping  A'TA  is  not  strongly  monotone  unless  the  matrix 
A* A  is  invertible.  Thus  convergence  of  iteration  (4)  is  not  guaranteed 
by  the  existing  theory  [1]  -  [3].  One  of  the  contributions  of  this  paper 
is  to  show  that  the  convergence  and  rate  of  convergence  properties  of 
Iteration  (4)  are  satisfactory  and  comparable  with  those  of  iteration  (8) . 
These  results  hinge  on  the  assumption  that  X  is  polyhedral,  and  it  is 
unclear  whether  and  in  what  form  they  hold  if  X  is  a  general  convex  set. 

In  Section  3  we  consider  a  classical  traffic  equilibrium  problem 
arising  in  several  contexts  including  communication  and  transportation 
networks,  which  can  be  modelled  in  terms  of  a  variational  Inequality  of 
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the  form  (1).  A  projection  algorithm  for  solving  this  problem  which  is 
essentially  of  the  form  (8)  has  been  given  by  Daferraos  [4].  Her  algorithm 
however  operates  in  the  space  of  link  flows,  and  involves  a  projection  iteration 
which  is  very  costly  for  large  networks.  We  consider  an  alternative  algorithm 
which  is  basically  of  the  form  (4)  and  operates  in  the  space  of  path  flows. 
Because  the  projection  iteration  can  be  carried  out  easily  in  this  space 
our  algorithm  is  much  more  efficient.  An  algorithm  which  has  several  similarities 
with  ours  has  been  proposed  by  Aashtiani  [5]  and  has  performed  well  in  compu¬ 
tational  experiments.  However,  Aashtiani' s  algorithm  cannot  be  shown  to 
converge  in  general.  By  contrast  the  results  of  Section  2  guarantee  convergence 
and  linear  rate  of  convergence  for  our  method.  There  are  also  other  methods 
[7] , [8] , [10] , [11] , [12]  for  solving  the  special  case  of  the  traffic  assignment 
problem  where  T  is  a  gradient  mapping  and  there  is  an  underlying  convex 
programming  problem.  Some  of  these  methods  [7],  [8],  [10]  are  of  the  projection^ 
type.  The  algorithm  of  the  present  paper,  however,  seems  to  be  the  first 
that  can  solve  the  general  problem,  is  suitable  for  large  networks,  and 


is  demonstrably  convergent. 


In  Section  4  we  consider  the  generalized  version  of  iteration  (4), 


Vl  '  PXk  t*k  -  “ 


where  {S^}  is  a  sequence  of  positive  definite  symmetric  matrices  with 
eigenvalues  bounded  above  and  bounded  away  for  zero.  Projection  algorithms 
of  this  type  include  Newton's  method  for  constrained  minimization  [14], [15],  and 
several  network  flow  algorithms  [5], [7]  -  [9].  Kxccpt  tor  the  case  where 
T  is  the  gradient  of  a  convex  function,  there  are  no  convergence  results 
in  the  literature  for  this  algorithm,  even  when  A  is  the  identity  matrix. 

We  show  that  if  care  is  taken  in  the  way  the  matrices  S^  are  allowed  to 


change,  Chen  Che  resulCing  algorithm  is  convergenc  at  a  rate  which  is  at 
least  linear.  We  also  provide  a  computational  example  Involving  a  traffic 
assignment  problem. 

2.  Projection  Methods  for  Variational  Inequalities 

* 

Let  X  be  the  set  of  all  solutions  of  the  variational  inequality 

(I)  .  We  have  that  X  is  polyhedral  and  is  given  by 

(9)  X*  -  {xeX  |  Ax  *=  y*} 

-ff 

where  y  is  the  unique  solution  of  the  variational  inequality  (6), 

For  any  xeRn  we  denote  by  p(x)  its  unique  projection  on  X*  with  respect  to 
the  norm  (5),  l.e. 

(10)  p(x)  -  arg  min  (||x-z||  |zeX*}  . 

We  recall  that  projection  on  a  convex  set  is  a  nonexpansive  mapping  (see 
e.g.  [13]),  so  we  have 

(II)  I  |p(xx)  -  P(x2>ll  £  "  X2H  »  V  • 

We  denote  the  mapping  A'TA  by  T,  i.e. 

(12)  T(x)  -  A’T(Ax)  ,  V  xeRn  . 

With  this  notation  we  have 

(13)  (x  -  x*)'  T(x*)  >  0,  V  xeX,  x*eX*  . 

We  also  denote  for  all  xeX,  a  >  0 

(14)  x(x,a)  -  P®  [x  -  as-1  T(x) ]  . 

With  this  notation  the  projection  iteration  (4)  is  written  as 

A 

(15)  -  x(xk,a),  xQeX. 

It  can  be  seen  that  we  can  also  obtain  x^+^  as  the  unique  solution  to  the  problem 

(16)  minimize  (*  -  x^) *  T(xfc)  +  |j(x  -  x^) '  S(x  -  xfc) 


subject  to  xeX 
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it  it 

In  view  of  (13)  it  is  easy  to  see  that  for  all  solutions  x  ex  and  a  >  0 
we  have 

(17)  x(x*,a)  -  x*  -  p(x*)  . 

* 

Thus  if  x^eX  for  some  k  the  algorithm  (15)  essentially  terminates. 

Our  main  result  is  given  in  the  following  proposition. 

Proposition  1:  Assume  that  T  is  LipschLtz  continuous  and  strongly  monotone, 
and  X  is  polyhedral. 

a)  There  exist  positive  scalars  q(S)  and  r(S)  depending  continuously  on  S, 
such  that  for  all  xeX  and  ct  >  0 

(18)  | |x(x,a)  -  p[x(x,a)J | |2  £  [1  -  2a  q(S)  +  a2r(S)]  ||x-p(x)|]2. 

b)  There  exists  a  >  0  such  that  for  all  ae(0,ci]  the  sequence  {x^} 
generated  by  Iteration  (15)  converges  to  a  solution  x*  of  the  variational 
inequality  (1).  The  rate  of  convergence'  is  at  least  linear  in  the  sense 
that  for  each  ae(0,a]  and  xq€X  there  exist  scalars  8  (depending  on  ct)and  q 

(depending  on  a  and  xq)  such  that  q  >  0,  Be (0,1),  and 

i  1^  “  x*l  i  i  S  3k  ,  k  -  0,1,... 

The  proof  of  Proposition  1  relies  on  the  following  lemma,  the  proof  of 
which  is  relegated  to  the  appendix.  The  lemma  is  easy  to  conjecture  in  terms 
of  geometrical  arguments  (see  Figure  1). 

Lemma  1:  Assume  the  conditions  of  Proposition  1  hold.  Let 

(19)  -X*  ■  (x£Rn| Ax  ■  y*) 

(20)  p(x)  -  arg  min  {j|x  -  zj  [  | zeX*}  . 

Then  p(x)  ■  p(p(x)]  and  there  exists  a  positive  scalar  0(S)  depending  continuously 
on  S  such  that  for  all  x£X 

(21)  |  |x  -  p(x)[  |  2  >_  n(S)|  j  p(x)  -  p(x)|  [  . 
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Proof  of  Proposition  1:  a)  Prom  (10)  we  obtain 

(22)  | |x(x,a)  -  p[Sc(x,a)]  1 |2  <_  [|x(x,ct)  -  p(x)  {  [ 2 
Also,  using  (17)  and  (11)  we  have 

(23)  ||£(x,a)  -  p(x)  i  | 2  -  |Jjt(x,a)  -  x[p(x)  ,at]  1 1 2 

<  III*  -  P(x)]  -aS_1[T(x)  -  T[p(x) ]]|  1 2  -  ||x  -  P(x)||2 

-  2a(T(x)  -  ftp(x)]]’  [x  -  p(x)]  +a2j  |s_1[T(x)  -  T[p(x)]]||2  . 

In  view  of  (22)  and  (23)  it  will  be  sufficient  to  show  that  there  exists 
positive  scalars  q(S)  and  r(S)  depending  continuously  on  S  such  that 

(24a)  [T(x)  -  f[p(x)]]'  [x  -  p(x)]  >  q(S)||x  -  p(x)||2  V  xcX 

(24b)  ■  |ls-1[T(x)  -  T(P(x)]]i|2  <  r(S)||x  -  p(x)||2  V  xCX  . 

To  prove  (24a)  we  first  use  (3)  and  (12)  to  obtain  for  all  xEX  and  y  =  Ax 

(25)  lT(x)  -  T[p(x)]]'(x  -  p(x) ]  »  [T(Ax)  -  T[Ap(x) ] ] ’ [Ax  -  Ap(x)] 

-  tT(y)  -  T(y*)]'(y  -  y*)  >  x  |y  -  y*|2  . 

We  now  consider  next  the  range  R(A)  and  the  nullspace  N(A)  of  A,  and  the 
orthogonal  complement 

N(A)+  -  { z6Rn |  z ' x  =  0,  V  xeN (A) }  . 

When  viewed  as  a  mapping  from  N(A)+  to  R(A),  A  is  one-to-one  and  onto,  and 
since  we  are  dealing  with  finite  dimensional  spaces  we  have  that  there 
exists  a  positive  scalar  0(S)  depending  continuously  on  S  such  that 

(26)  |w|2  >_  0(S)  | \z\ |2,  V  wER(A),  z£N(A)+,  w  -  Az. 

Since  for  all  y  m  Axf  and  p(x)  defined  by  (20),  we  have 

7  -  y*  -  A[x  -  p(x)] 

and  ^ 

x  -  p (x)  e  N(A)+ 


. 


we  obtain  from  (26) 
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(27)  |y  -  y*\2  >  0 (S) | |x  -  p(x)| |2. 

Using  the  Pythagorean  theorem  and  (21)  we  have 

(28)  | |x  -  p(x)| | 2  -  | | x  -  p(x)J | 2  +  | [p(x)  -  p(x)| | 2 

<  [i  +  ^sT  3  -  p(*)!I2- 

By  combining  (25),  (27),  and  (28)  we  obtain 

(T(x)  -  T[p(x)]]’[x  -  p(x)]  >  X8(S)  fl  ]_1|  |x  -  P(x)|l2 

so  (24a)  holds  with 

q(S)  -  X6(S)  [1  +  • 

To  show  (24b)  we  first  use  (2)  and  (12)  to  write 

(29)  ||s'1lT(x)  -  T[p(x) ]]j | 2  *  [T(y)  -  T(y*)]’  AS-1A'[T(y)  -  T(y*)] 

<  A(S)|T(y)  -  T(y*)|2  <  A(S)  L2|y  -  y*|2 

_i  » 

where  A(S)  is  the  largest  eigenvalue  of  AS  A  . 

By  the  same  argument  used  to  derive  (26)  we  can  assert  that  there  exists 
a  positive  scalar  y(S)  depending  continuously  on  S  such  that 

(30)  |w|2  <  Y(S) | 1 z | |2  ,  V  weR(A) ,  zeN(A)+,  w  -  Az  . 

Therefore 

(31)  |y  -  y*|2  £  y(S)!|x  -  P Cat)  1 1 2  £  Y(S)[|x  -  p(x)|i2  . 

Combination  of  (29)  and  (31)  yields  (24b)  for 

r(S)  -  A(S)L2Y(S). 

This  completes  the  proof  of  part  a) . 

b)  Let  a  be  any  positive  scalar  such  that 

a  < 


-29  (SI 

r(S) 


Then  using  (18)  we  have  for  all  ae(0,  a] 


|xk+l  "  p(xk+l)lf2  -  c«») I  Ix^  -  pC^II2 


where 


t(a)  -  1  -  2aq(S)  +  a  r(S)  . 

In  view  of  the  fact  a<_  3  <  ^gp'  it  is  easily  seen  that 


0  <  t(a)  <1  . 


We  have  from  (32) 


I  Ixjt  -  pCXfc)!  I2  i  t(a)k[  |xo  -  p(xo)[  |‘ 


and  using  (23)  and  (24)  we  obtain 


04)  ll^i-  P(*k)ll  <  tte)1/2||xk  -  pt^H  . 

By  the  triangle  inequality 

(35)  1  lxk+i  i’Hc+i  ’  1 1  +  I  l*k  '  p(xk}  1 1 


Combining  (35)  with  (33) ,  (34)  we  obtain 


\\+1  "  *kll  I  Cite'2  +  "  P<V 


<  [t(a)1/2  +  1J  t(a)k/2|[xo  -  p(xQ); 


8  -  t(a)‘ 


q  -  [t(a)  '  +  1]||*0  -  p(*0) 


Then  (36)  can  be  written  as 

1 1*^  -  \\  i  £  qBk  • 

For  all  k  >  0,  m  >  1  we  have 

(37>  iiv.  -  *ki  1 1 nv.  -  wiii+  ••  +iiv-i  •  *kH 

<5  <6l‘+m'1  +  . . .+  Bk>  - 

Since  in  view  of  t(a)  <  1  we  have  0  <  1,  it  follows  that  (x^)  is  a  Cauchy 

sequence  and  hence  converges  to  a  vector  x  .  Since  by  (33)  ix^  -  p(x^)} converges 

*  *  *  * 

to  zero  we  must  have  x  *  p(x  )  which  implies  x  GX  .  By  taking  limit  in 
(37)  as  m  -►  ®  we  obtain  for  all  k  *  0,1,... 

IK-x*ii  <  <»k 

where  q  •  1(1  -8)  .  Q.E.D. 

It  is  quite  remarkable  that  as  shown  in  Proposition  lb) ,  the 
sequence  { x^  }  converges  to  a  single  limit  point  even  though  the  set  of  solutions 
X  may  contain  an  infinite  number  of  points. 

3.  An  Algorithm  for  the  Traffic  Assignment  Problem 

In  this  section  we  consider  the  following  network  flow  problem. 

A  network  consisting  of  a  set  of  nodes  V  and  a  set  of  directed  links  L 
is  given,  together  with  a  set  W  of  node  pairs  referred  to  as  origin- 
destination  (0D)  pairs.  For  each  0D  pair  w€W  there  is  a  known  demand 
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dw>  0  representing  traffic  entering  the  network  at  the  origin  and  exiting 

at  the  destination.  For  each  OD  pair  w,  the  demand  d  is  to  be  distributed 

w 

among  a  given  collection  of  simple  directed  paths  joining  w.  We  denote  by 
Xp  the  flow  carried  by  path  p.  Thus  the  set  of  feasible  path  flow  vectors 
x  -  {Xp|pePw,  weW)  is  given  by 

(.38)  X  -  (x|  Z  x  «  d  ,  x  »  0,  V  pePw,  w€W)  . 

*P« 

Each  collection  of  path  flows  xeX  defines  a  collection  of  link  flows  y  , 
(i,j)eL  by  means  of  the  equation 

(39)  y  -  1  Z  5  (it j)  x  ,  V  (i,j) el 

weW  peP  p  p 

w 

where  6p(i,j)  ”1  if  path  p  contains  link  (i,j)  and  5p(i,j)  -  0  otherwise. 
The  vector  of  link  flows  y  *»  (y^j  |  (i»  j)eL)  corresponding  to  xeX  can  be 
written  as 

y  -  Ax 

where  A  is  the  arc-chain  matrix  defined  by  (39).  The  set  of  feasible  link 
flows  is  thus 

(40)  Y  -  AX  -  {y  1  y  -  Ax,  xeX> 

We  assume  that  for  each  link  (i,j)  L  there  is  given  a  function  T^  :  Y  -*•  R 
such  that  T^(y)  >  0  for  all  yeY.  The  value  of  T^  (y)  represents  a  measure 
of  delay  in  traversing  link  (i,j)  when  the  set  of  link  flows  is  y  (travel 
time  in  transportation  networks  [4],  [5],  marginal  delay  in  communication 
networks  [7],  [8]).  The  vector  with  components  T^ (y)  is  denoted  T(y). 

We  assume  that  T(y)  is  Llpschitz  continuous  and  strongly  monotone.  This 
is  a  reasonable  assumption  for  transportation  networks,  as  well  as  for 
communication  networks. 

For  each  xEX  and  corresponding  y  ■  Ax  the  vector  T(y)  defines  for 
each  w«W  and  Pe^w  a  length 
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(41)  T  (x)  -  Z  5  (i»j)  T .  ,  (y) 

P  (i.j)eL  P  iJ 

which  may  be  viewed  as  the  total  travel  time  of  path  p.  The  problem  is  to 
find  x*£X  such  that  for  all  Pe^w»  and  w£W 

(42)  T-(x*)  -  min  T  (x*)  ,  if  x-  >  0  . 

P  P  P 

pPw 

This  problem  is  based  on  the  user-optimization  principle  which  asserts 

that  a  traffic  network  equilibrium  is  established  when  no  user  may  decrease 

his  travel  time  by  making  a  unilateral  decision  to  change  his  route.  If  we 

denote  by  T(x)  the  vector  of  lengths  T  (x) ,  peP  ,  weW  then  it  is  easy  to 
'  p  w 

see  ([4],  [5])  using  (39)  and  (41)  that 

(43)  T(x)  -  A’T(Ax) 

and  that  the  problem  defined  earlier  via  (42)  is  equivalent  to  finding  a 

* 

solution  x  £X  of  the  variational  inequality 

(44)  (x  -  x*)'  T(x*)  ^  0,  V  xEX, 

which  is  of  the  form  (1) . 

Let  W  -  {w, ,w„, . . . ,wu}  and  consider  the  algorithm  of  the  previous 
l  z  n 

section  with  a  matrix  S  which  is  block  diagonal  of  the  form 


Here  each  matrix  corresponds  to  the  OD  pair  w^.  We  assume  that  each 
matrix  is  positive  definite  symmetric.  The  projection  iteration  (4) 
can  be  implemented  by  finding  solving  the  quadratic  program  (16). 

In  view  of  the  block  diagonal  form  of  S  and  the  decomposable  nature  of 
the  constraint  set  X  [ c . f .  (38) ]  the  quadratic  program  can  be  decomposed 


f 


into  a  collection  of  smaller  quadratic  programs  -  one  per  OD  pair  w€W. 

The  form  of  these  programs  for  the  case  where  each  matrix  is  diagonal 

with  elements  s  ,  p£P  along  the  diagonal  is 
P 

(45)  minimize  £  C(x  -  xS  T  (x.)  +  t£(x  -  x^)2} 

_ „  P  p  p  x  2a  p  p 


subject  to£  x-d  ,x>0,  V  p£P 
peP  p  wi  p  “  v 


where  x  pePw  i  *  1,...,M  are  the  components  of  the  vector  x.  .  Problem  (45) 

Involves  a  single  equality  constraint  and  can  be  solved  very  easily-essentially 

in  closed  form  [8],  [10].  The  convergence  and  linear  rate  of  convergence 

results  of  Proposition  1  apply  to  this  algorithm. 

The  preceding  algorithm  is  satisfactory  if  each  set  P  contains 

Wi 

relatively  few  paths.  In  some  problems  however  the  number  of  paths  in 

P  can  be  very  large  (for  example  P  may  contain  all  simple  paths  joining 
Wi  Wi 

w.)  .  In  this  case  it  is  preferable  to  start  with  a  subset  of  each  P  and 

augment  this  subset  as  necessary  as  suggested  by  Aashtiani  [5].  The 
corresponding  algorithm  is  as  follows : 

We  begin  with  a  subset  P°  C  P  for  each  w.EW  and  a  vector  of  initial 

Wi  wi  1 

path  flows  x  such  that  for  all  w.£W 
o  i 


x°  -  0 
P 


if  p*P 


At  the  kth  iteration  we  have  for  each  OD  pair  w  6W  a  subset  of  paths P  CP 

l  wt 

and  a  corresponding  vector  of  path  flows  x^  satisfying 

*p  -  0  if  P*V  • 


We  compute,  for  each  w.£W,  a  shortest  path  pK  £P  using  T. .(Ax.  )  as  length 

1  lj 
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for  each  link  (I,j).  We  set 


Bk+1 


P.  u(p..  > 


w. 


^  It 

(Note  here  that  p  may  already  belong  to  P  ) .  We  then  solve  the  quadratic 
1  wi 

programming  problem 


(46) 


nlol”lz'  ?fcfi  tc*P  -  =£  VV  +  iClp  -  xp)2} 


r-  -k 

If  x  , 
P 


(47) 


peP. 

subject  to  -S  x 

pePk+1  P 
w 

1 

,k+l 


dw±\  xpl°‘  ***< 


P 
k+1 


p£P  Is  the  solution  of  this  problem  we  set 

Wi  /  _k  k+1 

|xk  if  pCPw 

P  i 


k+1 


j  r  J-k+l 

lf  P^PW 

i 


Note  that  the  quadratic  programming  problem  (46)  is  the  same  as  (45) 

k+1 

except  for  the  fact  that  it  involves  only  paths  in  the  subset  P 

Wi 

The  subset  is  possibly  augmented  at  each  iteration  k  to  include 

lc 

the  current  shortest  path  p  .  The  expectation  here  is  that,  while  P 

Wi  Wi 

may  contain  a  very  large  number  of  paths,  the  actual  number  of  paths 

generated  and  Included  in  the  set  P  remains  small  as  k  increases. 

wi 

This  expectation  has  been  supported  by  computational  experiments  [5] . 

The  convergence  and  rate  of  convergence  properties  of  the  algorithm  (46) , 

(47)  are  identical  with  those  of  the  earlier  algorithm  based  on  (45)  as  the 

reader  can  easily  verify.  The  key  idea  is  based  on  the  fact  that  for  each 

k 

1,  P  contains  a  finite  number  of  paths  and  P  grows  monotonically  so 
W1  k  wi 

that  the  sequence  (P  }  converges  to  some  subset  of  paths  P  C  P 

wi  wi  W1 

By  applying  Proposition  1  it  follows  that  the  algorithm  converges  to  a 

solution  of  the  problem  obtained  when  P  is  replaced  by  P  .  Because 

wi  wi 
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P  is  augmented  at  each  iteration  with  the  current  shortest  path  it  is 
a  simple  matter  to  conclude  that  the  solution  x  obtained  satisfies  the 
required  minimum  travel  time  condition  (42). 

The  results  and  algorithms  of  this  section  can  be  strengthened 
considerably  in  the  case  where  T  is  in  addition  the  gradient  of  convex 
function  F  in  which  case  the  problem  is  equivalent  to  the  convex  programming 
problem  of  minimizing  F(y)  subject  to  ytY  or  minimizing  F(Ax)  subject  to  x€X. 
Under  these  circumstances  there  are  convergence  results  [8]  -  [10]  relating 
to  projection  type  algorithms  which  allow  for  a  variable  matrix  S  and 
for  coordinate  descent  type  iterations  whereby  each  iteration  is  performed 
with  respect  to  a  single  OD  pair  (or  a  small  group  of  OD  pairs)  and  all  OD 
pairs  are  taken  up  in  sequence.  In  fact  for  such  algorithms  it  seems  that 
it  is  easier  to  select  an  appropriate  value  for  the  stepsize  a.  Aashtiani's 
computational  experience  [5]  suggests  that  such  algorithms  also  work  well 
in  many  cases  where  T  is  not  a  gradient  mapping.  We  have  been  unable 
however  to  obtain  a  general  convergence  result  for  coordinate  descent 
versions  of  the  projection  algorithm.  The  possibility  of  changing  the 
matrix  S  from  one  iteration  to  the  next  is  considered  in  the  next  section. 

4.  An  Algorithm  with  Variable  Projection  Metric 

A  drawback  of  the  algorithms  of  Sections  2  and  3  is  that  the  matrix  S 
is  restricted  to  be  the  same  at  each  iteration.  Computational  experience 
with  optimization  problems  as  well  as  network  flow  problems  [5],  [6]  suggests 
that,  if  T  is  differentiable,  better  results  can  be  obtained  if  the 
matrix  S  is  varied  from  one  iteration  to  the  next  and  is  made  suitably 
dependent  on  first  derivatives  of  the  mapping  T  in  a  manner  which  approximates 


4 


Newton's  method.  We  have  not  been  able  to  show  the  result  of  Proposition  lb) 
for  algorithms  in  which  the  matrix  S  may  change  arbitrarily.  On  the  other 
hand  it  is  possible  to  construct  an  algorithmic  scheme  that  allows  for 
a  variable  matrix  S  but  at  the  same  time  incorporates  a  mechanism  that 
safeguards  against  divergence.  The  main  idea  in  this  scheme  is  to  allow 
a  change  in  the  matrix  S  only  when  the  algorithm  makes  satisfactory  progress 
towards  convergence.  The  algorithm  is  as  follows: 

A  set  S  of  positive  definite  symmetric  matrices  is  given.  It  is 

assumed  that  all  eigenvalues  of  all  matrices  ScS  lie  in  some  compact  Interval 

of  positive  real  numbers,  i.e.  there  exist  m^.n^  >  0  such  that 
2  «  2 

blJz|  _<  z  Sz  ^n^lzi  .for  all  z£R  >  SES.We  consider  the  algorithm 

<4fl>  Vi  ■  px“  l\  -  *<-k»  • 

where  S^eS  for  all  k  «  0,1,....  The  stepsize  a  is  such  that 


(49)  h  -  max  [1  -  2a  q(S)  +  Cl2  r(S)]  <  1, 

se S 

where  q(S)  and  r(S)  are  as  in  Proposition  la)  [c.f.  (18) J .  The  maximum  in  (49) 
is  attained  by  continuity  of  q(*)  and  r(»).  It  is  clear  that  there 
exists  ci  >  0  such  that  (49)  is  satisfied  for  all  ote(0,a].  Given  the 

matrix  Sfc+1  is  either  chosen  arbitrarily  from  S  or  else 

Sk+1  "  Sk 

depending  on  whether  the  quantity 

(50)  wk  -  (xk+1  -  V'  sk<Vn  -  V 

has  decreased  or  not  by  a  certain  factor  over  the  last  time  the  matrix  S  was 
changed.  More  specifically  a  scalar  Je(0,l) (typically  close  to  unity) 
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ls  chosen,  and  at  each  Iteration  k  a  scalar  is  computed  according  to 


gwk 

lf  uk  i  ”k 

"1*1  ■  j 

wk 

if  »k  >  ;k 

where  w^  is  given  by 

(50) 

and  initially  w^  ■  80 . 

We  select 

(52a)  Sk+1  - 

Sk 

*  lf  Wk+1 

wk  * 

(52b)  Vl*5 

•  if  *fcfl  < 

wk  * 

Thus  for  each  k,  the  scalar  represents  a  target  value  below  which 
w^  must  drop  in  order  for  a  change  in  S  to  be  allowed  in  the  next  iteration. 

We  first  show  that  if  {x^}  is  a  sequence  generated  by  the  algorithm 
Just  described  then 


(53)  lim  inf  w.  -  0 

k-*« 

Indeed  if  lim  inf  w.  >  0  then  S.  must  have  been  allowed  to  change  only 
lr*09  k  * 

a  finite  number  of  times  in  which  case  it  follows  from  Proposition  1  that 
.  .  * 

Ix^}  converges  to  a  solution  x  .  As  a  result  we  have  w^  0  contradicting 
the  earlier  assertion. 

Let  us  denote  by  I  1*1 1^  the  norm  corresponding  to  and  by  pk(z) 
the  projection  of  a  vector  z£Rn  on  X*  with  respect  to  j |*| |^.  We  have 
by  using  the  triangle  inequality,  (49),  (50)  and  Proposition  la) 

(54)  Wk  “I  !*k+l  “  Ik 

1  1 1\  “  pk°W  1  Ik  - 1  K+i  -  pk(W  1 1 
i  ll*k"  Pk^Hk”^  llxk-  pk(xk)!ik 

-  (1  -/h)  | -  Pk(*k) I ik 


Hence  (53)  Implies  lim  inf  | |x^  -  p^(x^) | j  ■  0  or  equivalently  (in  view 

•y  fj 

of  the  fact  >  0) 

(55)  U£*inf  l*k  "  pk(xk^  “  0  * 

This  means  that  at  least  a  subsequence  of  {x^}  converges  to  the  solution 

*  *  * 

set  X  .  We  will  show  that  in  fact  for  some  vector  x  EX,  and  some  scalars 

q  >0  and  £e(0,l)  we  have 

1^-  x*|  <  ,  6k  .  v  k  -  0,1.... 

.  * 

i.e.  tx.  }  converges  to  a  solution  x  at  a  rate  which  is  at  least  linear. 

We  have 

(56)  wk  -  I  Is*  S'V  +  P’/V  -  S+l1  Ik 

i  Us*  ■,k<xk)i|k+  iMv  -  vi'ik 
+  2|K*»k<*k>nk  nw-  viMk 

Using  (23),  (24),  and  (49)  we  obtain 

(57)  I  |pk<S>  -  Ski  I  Ik  i  h  I  K  *  S<xk>  I  Ik  • 

Combination  of  (56)  and  (57)  yields 

(58)  »k  <_  (1  +  i/h)2  I  -  PfcCKfc)  I  Ijf 

Also  in  view  of  (49)  and  Proposition  1  we  have  that  there  exists  a  scalar 
d  >_  1  such  that  for  all  k 

I  l*k+l  '  ’W*kH>  I  lk+1  -  I  l*k+l  *  Pkf*k+1}  I  lk+l 

id  Il*k+1  -  S'S+^Hk 

<  dh  iis  -  s<s>Hk 
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Inequalities  (63)  and  (65)  yield 


W2  , 

B»,.  <»,.,<  -  * 2- 1 -  ^  - 


kl  k2"1 


l-/h 


It  follows  that  if  k2  -  k^  1  then 


kj  -  <_  2  + 


g(l-  /h) 

<*(!+»£)  2 
In  h 


so  If  m  is  any  positive  integer  such  that 


we  have 


2  +  <5 

In  h  — 


k2  ”  ^i  £® 


for  any  two  successive  indices  k^  k^  in  K  .  It  follows  using  (62) 


.  /S  m  .k 

w,  <  w  (0  ) 

k  —  o 


V  kcK  . 


Using  (65),  (68)  and  (69)  we  also  obtain 


d(l+i/h)2  w  -  .  - 

v  <  - =— 2.  (g  m  )k-“ 

wk  —  1  VC  VP  ' 


V  k  t< 


Combining  (69)  and  (70)  we  have  for  some  scalar  q  >  0 

1 


wv  i,  °  )l 


¥  k  -  0,1,, 


(71) 
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Since 

■1^1“  *ki2i  “wk 

relation  (71)  yields 

1 

Since  0  <  1,  it  follows  in  exactly  the  same  manner  as  in  the  proof  of 
Proposition  lb)  that  {x^}  a  Cauchy  sequence  which  converges  to  a 
vector  x  6  X  .  Furthermore  for  some  q  >0  and  0e(O,l)  we  have 

(72)  \x^  -  x*j  <  q0k  ,  k  -  0,1,... 

We  have  thus  proved  the  following  proposition: 

Proposition  2:  There  exists  a  >  0  such  that  if  ct6(0,a],  a  sequence 
{x^J  generated  by  iteration  (48)  with  {S^}  selected  according  to  (50)  -  (52) 
converges  to  a  solution  x*  of  the  variational  inequality  (1) .  Furthermore 
there  exist  scalars  q  >  0  and  0e(_O,l)  such  that 

1^  -  x*|  <  qek  ,  * 

5.  A  Computational  Example 

In  this  section  we  report  computational  results  for  a  traffic 
assignment  problem.  The  corresponding  network  is  shown  in  Figure  2,  and 
may  be  viewed  as  a  model  of  a  circular  highway.  There  are  five  origins 
and  destinations  numbered  1, 2,3,4, 5  and  connected  through  the  highway 
via  entrance  and  exit  ramps.  We  consider  the  five  0D  pairs  (1,4),  (2,5), 
(3,1),  (4,2),  (5,3).  Each  0D  pair  has  two  paths  associated  with  it  - 
the  clockwise  and  counterclockwise  paths  on  the  corresponding  circle. 

The  expressions  for  the  travel  time  on  each  link  are  shown  in  Figure  2 


1.  Highway  Links:  17,27,37,47,57,18,28,38,48,58 

2.  Exit  Ramps:  14,24,34,44,54,12,22,32,42,52 

3.  Entrance  Ramps:  11,21,31,41,51,13,23,33,43,53 

4.  Bypass  Links:  15,25,35,45,55,16,26,36,46,56 
Delay  on  Links  [where  g  is  defined  by  g(x)  *  1  +  : 

1.  Delay  on  Highway  Link  k:10.g[Flow  on  k]  1-  2*yg[Flow  on  exit  ram  from  k] 

2.  Delay  on  Exit  Ramp  k:  g[Flow  on  k] 

3.  Delay  on  Entrance  Ramp  k:g  [Flow  on  k]  +  y g  [Flow  on  Bypass  Link 

merging  with  k] 

4.  Delay  on  Bypass  Link  k:  g(Flow  on  k] . 

Remark:  Flow  is  not  allowed  to  use  exit  ramp  not  leading  to  its  destination 


where  the  function  g  is  given  by 


g  00  -  1  +  x  +  x2 

Different  values  of  the  nonnegative  scalar  y  represent  different  degrees  of 
dependence  of  the  travel  times  of  some  links  on  the  flows  of  other  links. 
The  problem  is  equivalent  to  an  optimization  problem  if  and  only  if  there 
is  no  such  dependence  (y=0) . 

Tables  1  and  2  list  representative  computational  results  for  two 

✓ 

demand  patterns,  three  values  of  y,  and  fifteen  iterations  of  two  different 
algorithms  labeled  "all-at-once"  and  "one-at-a-time"  and  described  below. 
The  number  shown  in  the  tables  for  each  iteration  k  is  the  following 
normalized  measure  of  convergence 


Ax 


(73) 


w,k 


AT 


w,k 


all  0D  pairs  w  w 


min,w,k 


where  d  is  the  demand  of  the  OD  pair  w,  Ax  .is  the  portion  of  the  demand 

W  W)K 

that  does  not  lie  on  the  shortest  path  of  the  OD  pair  w  at  the  end 

of.  iteration  k,  AT  ,  is  the  difference  of  the  travel  times  of  the  longest 

W 1 1C 

and  shortest  paths  and  T  .  ,  is  the  travel  time  of  the  shortest  path. 

Clearly  the  expression  (73)  is  zero  if  and  only  if  the  corresponding  traffic 

assignment  is  optimal.  The  starting  flow  pattern  in  all  runs  was  the  worst 
possible  whereby  all  the  demand  of  each  OD  pair  is  routed  on  the  counter¬ 


clockwise  path.  The  results  suggest  that  the  algorithms  yield  near  optimal 
flow  patterns  after  very  few  iterations  and  subsequently  continue  their 
progress  at  a  fairly  satisfactory  rate.  This  type  of  convergence  behavior 
is  consistent  with  the  one  observed  for  related  algorithms  tested  in  [6]. 

The  "all-at-once"  algorithm  is  the  one  of  the  previous  section 
[c.f. (48),  (50)  -  (52) J  with  the  projection  matrix  being  diagonal 


TABLE  Is  Demands  d(l,4)  -  0.1;  d(2,5)  -  0.2;  d(3,l)  -  0.3;  d(4,2)  -  0.4;  d(5,3)  -  0.5 


TABLE  2:  Demands  d(l,4)  *  1,  d(2,5)  *  8,  d(3,l)  *  1,  d(4,2)  ■  8,  4(5,3) 
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1 

[c.f.  (45)].  For  each  iteration  k  for  which  S.  was  allowed  to  change 

K  ! 
[c.f. (50)  -  (52)]  the  diagonal  element  of  S,  corresponding  to  any  one  t 

k  '  t 

path  was  taken  equal  to  the  sum  of  the  first  derivatives  of  the  travel  1 

r 

times  of  links  on  that  path  evaluated  at  the  kth  flow  x,  .  This  corresponds 

*  *  j 

to  a  diagonal  approximation  of  Newton's  method  (c.f. [8],  [9]).  As  a  result  I 

stepslzes  near  unity  typically  give  satisfactory  convergence  behavior.  In  J 

( 

all  runs  we  used  a  stepsize  a  -  0.8  which  is  probably  a  bit  on  the  high  j 

side.  The  scalar  5  used  in  the  test  for  allowing  the  matrix  to 

change  [c.f. (51)]  was  taken  to  be  0.99  in  all  runs.  One  of  the  most 

■ 

interesting  observations  from  our  experimentation  was  that  this  test  was 
passed  at  every  iteration  and  so  the  matrix  was  changed  at  every  iteration. 

The  "one-at-a-time"  algorithm  is  similar  to  the  "all-at-once"  algorithm 
described  above.  The  only  difference  is  that  the  projection  is  carried 
out  with  respect  to  a  single  OD  pair  with  flows  corresponding  to  the  other  j 

OD  pairs  being  kept  fixed,  and  all  OD  pairs  are  taken  up  in  sequence. 

An  iteration  consists  of  a  cycle  of  five  projection  subiterations  (one 
per  OD  pair).  Each  subiteration  is,  of  course,  followed  by  reevaluation 
of  the  travel  time  of  each  link.  Algorithms  of  this  type  resemble  coordinate 
descent  methods  for  unconstrained  optimization  and  have  been  suggested  in 
the  context  of  network  flows  in  [5]  and  [8]— [10] .  The  stepsize  a  was  taken 
to  be  unity  in  all  runs.  Also  $  was  chosen  to  be  0.99  and  again  it  turned 
out  that  the  matrix  was  allowed  to  change  in  every  iteration.  It  is 
interesting  to  note  that  for  this  stepsize  and  this  particular  example 
the  algorithm  tested  is  equivalent  to  an  algorithm  in  the  class  proposed  by 
Aashtiani  [5].  Tables  1  and  2  indicate  a  better  performance  for  the 
"one-at-a-time”  algorithm.  However,  there  is  no  convergence  proof  available 
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for  this  algorithm  at  present  except  In  the  case  where  the  corresponding 
variational  Inequality  is  equivalent  to  a  convex  optimization  problem. 

Appendix  :  Proof  of  Lemma  1 

Ue  have  by  the  Pythagorean  theorem 

II*  -  p(*)ll2  +1  |p(*>  -  p(*)II2  -  II*  -  p<*>ll2l  II*  -  ptpOOJ  1 12 

-  1 1*  -  p(*)| 1 2  +  I |p(*>  -  p[p(*)J II2 

Hence 

(A.  1)  ||p(x)  -  p (x)  |  | 2  <_  |  |p00  -  p[pCx)]||2. 

Since  p[p(x)]  by  definition  is  the  unique  solution  of  the  problem  of 
minimizing  ||p(x)  -  z||2  over  xcX  ,  and  p(x)eX  ,  it  follows  that  equality 
holds  in  (A.l)  and  p(x)  -  p[p(x)]. 

Let  p^(x)  and  p^(x)  be  the  projections  of  x  on  X  and  X  respectively 
relative  to  the  standard  norm  [*|  ,  i.e.  p^(x)  ■  p(x),  p^(x)  ■  p(x)  for 
S  equal  to  the  identity.  In  order  to  show  the  existence  of  a  continuous 
T)(S)  such  that  (21)  holds  it  will  suffice  to  show  the  existence  of  a  scalar 
nx  >  0  such  that  for  all  x£X 

(A.2)  |x  -  ^(x)  |2  >,  Tlj^lpj^Cx)  -  px(x)  [2  . 

In  order  to  see  this  let  A(S)  and  X(S)  be  the  largest  and  smallest 
eigenvalue  of  S.  We  have  for  all  x,  zcRn 

X  (S)  |x  -  z\2  <_  ||x-  z|  |2  <  A(S)  |  x  —  z|2  . 
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It  follows  Chac 

<A,3)  Ycsyi  1*  *  p(x) 1 12  i  ix "  Pi^x) I2 

(a.4)  acsH^*  ’  p(*)ll2  l  U  -  pi(*)l2  • 

By  the  Pythagorean  theorem  we  have 

Ipx(x)  -  p^x)]2  -  [x  -  px(x)|2  -  |x  -  ?x(x)|2 
so  (A. 2)  can  ve  rewritten  as 

(A-5)  |x  -  px(x)|2  >  |x  -  p1(x)|2. 

From  (A. 3),  (A. 4),  and  (A. 5)  we  obtain 

(A.  6)  |  |x  -  p(x)  |  | 2  >_  Y(S)  |  |x  —  p(x)  |  j 2 


where 


Y«)  ■  ^ 


l-hlj  A(s) 


By  the  Pythagorean  theorem  we  have 

|  |x  -  p(x)  ||2  -|  |x  -  p(x)  I  I 2  +  I  |p(x)  -  p(x)  I  1 2 
and  by  using  this  equation  In  (A. 6)  we  obtain  the  desired  relation 
|  |x  -  p(x)  | \2  >  n(S)  I  |p(x)  -  p(x)  j  |2 


with 


n(s)  - 


rtf) 


1  -  Y(S) 


We  now  show  existence  of  an  >  0  such  that  (A. 2)  holds  for 
all  xeX.  By  the  preceding  analysis  this  Is  sufficient  to  prove  the  lemma. 

For  each  xeX  consider  the  tangent  cone  Cx  of  X  at  p(x),  i.e.  the  set 


(A.  7) 


Cx  ■  {*|  there  exists  a  >  0  such  that  Ip 00  +  azjcx}. 
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Let  #  be  Che  collection 

(A.  8)  -  (cj  xeX}  . 

It  is  easily  seen  that  ^  is  a  finite  collection,  i.e.  for  some  finite 
set  J  C  X  we  have 

(A.  9)  {Cj|jej} 

Indeed  since  X  is  polyhedral  it  can  be  represented  by  definition  ([16], 

Section  19),  in  terms  of  a  finite  number  of  vectors  d.,...,d  cRr'  and 

X  m 

scalars  bt . b  as 

1  m 

X  ■  {xld^  x  £  bif  »m) 


It  is  easily  seen  that 

r  i  '  i 

cx  "  i*l  djZ  £  0,  V  i  -  1 . m  such  that  dip(x)  *  b^}  . 

Clearly  there  is  only  a  finite  number  of  sets  of  the  above  form.  In 
what  follows  Pjj(z)  denotes  the  projection  of  a  vector  zeRn  on  a  closed 
convex  set  flCRn  with  respect  to  the  standard  norm  |»|  .  The  essence  of 
the  proof  of  Lemma  1  is  contained  in  the  following  lemma. 


Lemma  A.l:  For  jeJ  let 

(A.  10)  M*  -CjHn(A) 


(A.  11) 


Z 

j 


-  lz\  zee,,  P  *(z)  -  0}  , 

j  Mj 


where  N(A)  is  the  nullspace  of  A.  Then  for  each  JeJ  there  exists  a  scalar 
>  0  such  that 

<A*12>  I*  “  P»CA)**^2  -  njl*l2  V  seZj  . 
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Proof ;  Assume  Che  contrary,  i.e.  chat  there  exist  jeJ  and  sequences 
Cz^CZj  »  (Hj)C  R  such  that 

<*•«>  l*k*  W*k)|2<nJ|,k|2  •  "i*0- 

We  then  have  z^  i*  0,  V  k,  and  since  both  Mj  and  Zj  are  clearly  cones  with 
vertex  at  the  origin  we  can  assume  that  |z^[  -  1,  V  k.  Let  z  be  a  limit 
point  of  {z^}  .  We  have  |zj  “1  and  by  taking  limit  in  (A*  13)  we  obtain 


i.e.  zeN (A) . 
follows  that 


Since  Cj  is  closed  we  also  have  zeCj  and  hence  zeM* 


It 


(A.  14) 


V  <*> 

Mj 


On  the  other  hand  since  z.eZ,  we  have  P„*(z.  )  *  0,  V  k  which  implies 

x  J  Mj  K 

that  P  *(z)  ■  0.  Since  |z|  ■  1,  this  contradicts  (A. 14).  Q.E.D. 

We  now  show  that  the  desired  relation  (A. 2)  holds  with 

(A.  15)  -  min  {rijljeJ}  >  0  . 

Choose  any  xeX  and  let  jeJ  be  such  that 

cx  -  Cj  . 

Let 

(A.  16)  z  »  x  -  Pj^Cx)  . 

By  a  simple  translation  argument,  the  fact  that  P^(x)  is  the  projection  of  x 
on  x*H  {x|Ax»y*} implies  that  the  projection  of  z  on  M*  is  the  origin  so 

that  zCZj.  A  similar  argument  shows  that 
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(A.17)  Pj^C*)  -  P^(x)  ■  Pn(a)^  • 

Using  (A.  16)  and  (A. 17)  in  (A. 12)  we  obtain 

|x  -  p£x)|2  >  Hj|x  -  PL(x) | 2 

and  (A. 2)  follows  from  the  definition  (A. IS)  of  Q.E.D. 
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